\documentclass{article}
\usepackage{geometry}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

% useful packages.
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{enumerate}
\usepackage{graphicx}
\usepackage{multicol}
\usepackage{fancyhdr}
\usepackage{layout}
% \usepackage{ctex}
\usepackage{listings}
\usepackage{subfigure}
\usepackage{setspace}

% some common command
\newcommand{\dif}{\mathrm{d}}
\newcommand{\avg}[1]{\left\langle #1 \right\rangle}
\newcommand{\difFrac}[2]{\frac{\dif #1}{\dif #2}}
\newcommand{\pdfFrac}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\OFL}{\mathrm{OFL}}
\newcommand{\UFL}{\mathrm{UFL}}
\newcommand{\fl}{\mathrm{fl}}
\newcommand{\op}{\odot}
\newcommand{\Eabs}{E_{\mathrm{abs}}}
\newcommand{\Erel}{E_{\mathrm{rel}}}
\newcommand{\RNum}[1]{\uppercase\expandafter{\romannumeral #1\relax}}

\usepackage{xcolor}
\usepackage{fontspec} 
\definecolor{dkgreen}{rgb}{0,0.6,0}
\definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{comment}{rgb}{0.56,0.64,0.68}

\newfontfamily\monaco{Monaco}
\lstset {
aboveskip=3mm,
belowskip=3mm,
showstringspaces=false,       % underline spaces within strings
columns=flexible,
framerule=1pt,
rulecolor=\color{gray!35},
backgroundcolor=\color{gray!5},
basicstyle={\small\monaco},           % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\tiny\monaco\color{gray},  % the style that is used for the line-numbers
numbersep=5pt,                  % how far the line-numbers are from the code
commentstyle=\color{comment},
keywordstyle=\color{blue},
stringstyle=\color{dkgreen},
tabsize=2,                      % sets default tabsize to 2 spaces
captionpos=b,                   % sets the caption-position to bottom
breaklines=true,                % sets automatic line breaking
breakatwhitespace=false,        % sets if automatic breaks should only happen at whitespace
escapeinside={\%*}{*)},            % if add LaTeX within your code
morekeywords={*,...}               % if add more keywords to the set
}


\begin{document}
\title{Homework \#6}
\pagestyle{fancy}
\lhead{Name Li HuiTeng 3180102114}
\chead{ NumAnalysis\#6}
\rhead{Date 21.12.14}

\section{Theoretical questions}

\RNum{1}
\begin{proof}
For (a), denote $y_0=y(-1),y_1=y(0),y_2=y(1),y_1'=y'(0)$. Apply Hermite interpolation and 
\begin{align*}
  \begin{array}{c|cccc}-1&y_0&&&\\0&y_1&y_1-y_0&&\\0&y_1&y_1'&y_1'-y_1+y_0&\\1&y_2&y_2-y_1&y_2-y_1-y_1'&\frac{y_2-y_0-2y_1'}2\end{array}\Rightarrow \\
  p_3(x)=y_0+(y_1-y_0)(x+1)+(y_1'-y_1+y_0)x(x+1)+\frac{y_2-y_0-2y_1'}2x^2(x+1).
\end{align*}
By integration, 
\begin{align*}
  \int_{-1}^1p_3(x)dx&=2\int_0^1y_0+y_1-y_0+(y_1'-y_1+y_0+\frac{y_2-y_0-2y_1'}2)x^2dx\\
  &=2(y_1+\frac{y_2+y_0-y_1}6)\\
  &=\frac{1}{3}(f(-1)+4f(0)+f(1)),
\end{align*}
which yields $\textbf{$\omega$}_0=\{\frac{1}{3},\frac{4}{3},\frac{1}{3}\}$.

For (b), by Theorem 3.33, 
\begin{align*}
  E^S(y)&=\int_{-1}^{1}R_3(y;x)dx\\
  &=\int_{-1}^{1}\frac{f^{(4)}(\xi(x))}{24}(x-1)x^2(x+1)dx\\
  &=-\frac{1}{90}f^{(4)}(\zeta ),
\end{align*}
where $\xi(x),\zeta \in [-1,1]$ and the third step follows from integral mean value theorem due to the positivity of $(x-1)(x+1)x^2$.

Now notice when $f \in \mathbb{P}_3[x]$, $E^S(y)=0$. Therefore, such rule is accurate for $\mathbb{P}_3[x]$.

For (c), denote $t=\frac{2}{b-a}(x-a)-1$, and we have 
\begin{align*}
  \forall p(x) \in \mathbb{P}_3[x],\\
  \int_{a}^{b}p(x)dx&=\int_{-1}^{1}p(\frac{b-a}{2}t+\frac{b+a}{2})\frac{b-a}{2}dt\\
  &=\frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b)),
\end{align*}
where the first step follows from a change of variable and the second from (a)'s result and (b)'s fact that such rule is accurate for all cubic polynomials. This yields $\omega=\{\frac{b-a}{6},2\frac{b-a}{3},\frac{b-a}{6}\}$.

For error estimation, we have 
\begin{align*}
  E^S(y)&= \int_a^bR_3(x)dx\\
&=\int_{a}^b\frac{f^{(4)}(\xi(x))}{24}(x-a)(x-\frac{a+b}{2})^2(x-b)dx\\
&=(\frac{b-a}{2})^5\int_{-1}^{1}\frac{f^{(4)}(\xi(x(t)))}{24}(t-1)t^2(t+1)dt\\
&=-\frac{(b-a)^5}{2880}f^{(4)}(\zeta),
\end{align*}
where $\xi(x),\zeta \in [a,b]$, the third step follows from a change of variable and the fourth step is just same as (b)'s third step. 
\end{proof}

\RNum{2}
\begin{proof}
Since
\begin{align*}
 f'&= -2*x*exp(-x^2)\\
 f''&= 4*x^2*exp(-x^2) - 2*exp(-x^2)\\
 f'''&= 12*x*exp(-x^2) - 8*x^3*exp(-x^2)\\
  f''''&=12*exp(-x^2) - 48*x^2*exp(-x^2) + 16*x^4*exp(-x^2),
\end{align*} 
then plot $f''$ and $f''''$. 

\includegraphics[width=0.45\textwidth]{2_df2.png}
\includegraphics[width=0.45\textwidth]{2_df4.png}

Therefore it's clear that $\|f''\|_{\infty}=|f''(0)|=2$ and $\|f''''\|_{\infty}=|f''''(0)|=12$.

Also, since $0<\int_0^1exp(-x^2)dx<\int_0^1exp(-x)<1$, in order to preserve 6 correct decimal places, 
the absolute error should be no greater than $0.5\times 10^{-6}$.

According to Theorem 6.18, for composite trapezoidal rule, we have 
$$|E_n^T| \leq \frac{1}{12}\|f''\|_{\infty}h^2=\frac{h^2}{6} .$$
Hence there should be 
$$\frac{h^2}{6} \leq 0.5\times 10^{-6}.$$
By calculation, we have $n=\frac{1}{h}\geq 577.35027$, and $N^T=578$.

According to Theorem 6.20, for composite Simpson's rule, we have 
$$|E_n^S| \leq \frac{1}{180}\|f''''\|_{\infty}h^4=\frac{h^4}{15} .$$
Hence there should be 
$$\frac{h^4}{15} \leq 0.5\times 10^{-6}.$$
By calculation, we have $n=\frac{1}{h}\geq 19.10886$, and $N^S=20$.
\end{proof}

\RNum{3}
\begin{proof}
For (a), it's equivalent to $(1,\pi)=(x,\pi)=0$. These two conditions yields 
\begin{align*}
  &\int_{0}^{\infty}(t^2+at+b)e^{-t}=2+a+b=0,\\
  &\int_0^{\infty}(t^2+at+b)te^{-t}=6+2a+b=0.
\end{align*}
Hence $a=-4,b=2$. And the orthogonal polynomial is 
$\phi_2(t)=t^2-4t+2$ with its zeros at $t_1=2+\sqrt{2},t_2=2-\sqrt{2}$.\\
For (b), since Gauss quadrature is exact for all constants and linear polynomials, 
we have 
\begin{align*}
  \omega_1+\omega_2 &= \int_0^{\infty}e^{-t}dt=1,\\
  t_1\omega_1+t_2\omega_2 &= \int_0^{\infty}te^{-t}dt=1,
\end{align*}
which yields $\omega_1=\frac{2-\sqrt2}{4},\omega_2=\frac{2+\sqrt2}{4}$.

By Theorem 6.36, 
\begin{align*}
  E_2(f)&=\frac{f^{(4)}(\tau )}{24}\int_0^{\infty}\pi_2(t)^2e^{-t}dt=\frac{f^{(4)}(\tau)}{6},
\end{align*}
where $\tau > 0$.

For (c), apply formula (b) and 
\begin{align*}
  I_{GL}&=\omega_1\frac{1}{1+t_1}+\omega_2\frac{1}{1+t_2}\\
&=\frac{4}{7},\\
|E_2(f)|&\leq\frac{\|f^{(4)}\|}{6}=1, 
\end{align*}
where the last step follows from $\|f^{(4)}\|_{\infty}=\|6(t+1)^{-4}\|_{\infty}=6$.

The true absolute error is 
\begin{align*}
  E_{abs}=|I_{GL}-I|=2.49188\times10^{-2},
\end{align*}
which is much smaller than the estimation bound. Hence $(\tau+1)^{-4}=E_{abs}$ yields $\tau=1.51691$.
\end{proof}

\RNum{4}
\begin{proof}
For (a), since 
\begin{align*}
  p(x_k)&=\sum_{m=1}^n[h_m(t)f_m+q_m(t)f'_m]=h_k(x_k)f_m+q_k(x_k)f'_k=(a_k+b_kx_k)f_k+(c_k+d_kx_k)f_k'\\
  p'(x_k)&=\sum_{m=1}^n[b_ml_m^2+2(a_m+b_mt)l_ml_m']f_m+[d_ml_m^2+2(c_m+d_mt)l_ml_m']f'_m\\
  &=[b_k+2(a_k+b_kx_k)l_k'(x_k)]f_m+[d_k+2(c_k+d_kx_k)l_k'(x_k)]f_m' 
\end{align*}
we have $\forall m =1,\cdots,n$,
\begin{align*}
  \left\{\begin{array}{l}a_m+b_mx_m=1\\c_m+d_mx_m=0\\b_m+2(a_m+b_mx_m)l_m'(x_m)=0\\d_m+2(c_m+d_mx_m)l_m'(x_m)=1\end{array}\right.\Rightarrow\left\{\begin{array}{l}a_m=1+2l_m'(x_m)x_m\\b_m=-2l_m'(x_m)\\c_m=-x_m\\d_m=1\end{array}\right.
\end{align*}
For (b), integrating the above polynomial over the interval $[a,b]$ for the weight function $\rho$ yields 
\begin{align}
  \int_a^b p(x)\rho(x)dx = \sum_{k=1}^{n}(\omega_kf_k+\mu_k f_k'), 
\end{align}
where $\omega_k$ is dependent on $a_k,b_k,\rho$ and the interval, $\mu_k$ on $c_k,d_k,\rho$ and the interval, i.e.
\begin{align*}
  \left\{\begin{array}{l}\omega_m=\int_a^b(a_m+b_mt)l_m(t)^2\rho(t)dt\\\mu_m=\int_a^b(c_m+d_mt)l_m(t)^2\rho(t)dt\end{array}\right. .  
\end{align*}

Next we show such formula is accurate for all $f(x) \in \mathbb{P}_{2n-1}$ by calculating $E_n(f)$.

$\forall f \in \mathcal{C}^{2n}$, apply (a)'s Hermite
interpolation. Then the error term of Hermite interpolation is 
$$\forall x \in [a,b],f(x)-p(x)=\frac{f^{(2n)}(\xi(x))}{(2n)!}\prod_{k=1}^n(x-x_k)^2,\xi \in [a,b]$$ 
After integration, by mean value theorem of integration, 
$$E_n(f)=\int_a^b\frac{f^{(2n)}(\xi(x))}{(2n)!}\prod_{k=1}^n(x-x_k)^2=\frac{f^{(2n)}(\tau)}{(2n)!}\int_a^b\prod_{k=1}^n(x-x_k)^2,\tau \in [a,b].$$

Hence for all $f(x) \in \mathbb{P}_{2n-1}$, $f^{(2n)}=0$ and $E_n(f)=0$. 

For (c), if one requires that after integration all coefficients of the derivative
terms should be zero, then we have 
$$
\forall m=1,\cdots,n, \mu_m=\int_a^b(t-x_m)l_m(t)^2\rho(t)dt=0
$$
\end{proof}

\RNum{5}
\begin{proof}

  For (a), expand $u(x+h),u(x-h)$ at x:
	\begin{align*}
		u(x+h) & =u(x)+hu'(x)+\frac{h^2}{2}u''(x)+\frac{h^3}{6}u'''(x)
		+\frac{h^4}{24}u^{(4)}(\xi_1),\xi_1 \in [x,x+h],                              \\
		u(x-h) & =u(x)-hu'(x)+\frac{h^2}{2}u''(x)-\frac{h^3}{6}u'''(x)
		+\frac{h^4}{24}u^{(4)}(\xi_2),\xi_2 \in [x-h,x],\\\Rightarrow
    D^2(x)&=\frac{u(x+h)+u(x-h)-2u(x)}{h^2}=u''(x)+\frac{h^2}{12}u^{(4)}(\mu),\mu \in [\xi_2,\xi_1],
  \end{align*}
  where the last step follows from mean value theorem of continuous function. 

  Since $u \in \mathcal{C}^4$, such method is second order accurate.

  If one consider perturbed errors $\epsilon \in [-E,E]$, the formula turns to 
  \begin{align*}
    |\widetilde{D}^2(x)-u''(x)|&=|\frac{\widetilde u(x+h)+\widetilde u(x-h)-2\widetilde u(x)}{h^2}|\\
    &\leq \frac{|u(x+h)+u(x-h)-2u(x)|+4E}{h^2}=\frac{h^2}{12}|u^{(4)}(\mu)|+\frac{4E}{h^2},\mu \in [x-h,x+h].
  \end{align*}
  By Inequality of arithmetic and geometric means, the error bound $E$ satisfies 
  \begin{align*}
    E\geq 2\sqrt{\frac{h^2}{12}|u^{(4)}(\mu)|\frac{4E}{h^2}}=2\sqrt{\frac{|u^{(4)}(\mu)|E}{3}},
  \end{align*}
  where the equality holds if and only if $h=\sqrt[4]{\frac{48E}{|u^{(4)}(\mu)|}}$, supposed the fourth order derivative is nonzero.

  For (b),  
  \begin{align*}
    U_0&=2u(x),\\
    U_1=u(x+h)+u(x-h)&=2[u(x)+\frac{h^2}{2}u''(x)+\frac{h^4}{24}u^{(4)}(x)+\frac{h^6}{720}u^{(6)}(x)]+O(h^8),\\
    U_2=u(x+2h)+u(x-2h)&=2[u(x)+4\frac{h^2}{2}u''(x)+16\frac{h^4}{24}u^{(4)}(x)+64\frac{h^6}{720}u^{(6)}(x)]+O(h^8),
  \end{align*}
  to solve out $a,b,c$ such that $\frac{aU_1+bU_2+cU_0}{h^2}-u''(x)=O(h^4)$, we have 
  \begin{align*}
    \left\{\begin{array}{l}2(a+b+c)=0\\2(a/2+4b/2)=1\\2(a/24+16b/24)=0\end{array}\right.\Rightarrow\left\{\begin{array}{l}a=\frac43\\b=-\frac1{12}\\c=-\frac54\end{array}\right.  .  
  \end{align*} 

  Hence we obtain a five point stencil with accuracy 4 such that 
  \begin{align*}
    D^2_4(x)=\frac{1}{12}[-u(x-2h)+16u(x-h)+15u(x)+16u(x+h)-u(x+2h)].
  \end{align*}

  Similarly, we have 
  \begin{align*}
    D^2_4(x)=u''(x)-\frac{1}{90}h^4u^{(6)}(\mu),\mu \in [x-2h,x+2h]. 
  \end{align*}

  Again assuming perturbed errors $\epsilon \in [-E,E]$, the formula turns to 
  \begin{align*}
    |\widetilde{D}^2_4(x)-u''(x)|\leq \frac{1}{90}h^4u^{(6)}(\mu)+\frac{49E}{12h^2}.
  \end{align*}
  Denote $|u^{(6)}|=U_6$, and the error bound $E_4$ satisfies 
  \begin{align*}
    E_4=\frac{1}{90}h^4U_6+\frac{49E}{24h^2}+\frac{49E}{24h^2}\geq 3(\frac{1}{90}h^4U_6\frac{49E}{24h^2}\frac{49E}{24h^2})^{1/3}=3(\frac{2401U_6E^2}{51840})^{1/3}
  \end{align*}
  where the equality holds if and only if $\frac{1}{90}h^4U_6=\frac{49E}{24h^2}$, i.e. $h=(\frac{735E}{4U_6})^{1/6}$ .
\end{proof}

\end{document}
%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
